1. Field of the Invention
The present invention generally relates to the modeling of hydrocarbon reservoirs and particularly to providing a more consistent and accurate representation of fluid flows and pressure gradients.
2. Description of the Related Art
Locating petroleum deposits, drilling wells and managing the extraction of petroleum from a deposit is done at great expense. Due to the complicated and critical nature of locating and managing hydrocarbon deposits, new professions have evolved and special technologies have been developed in order to reduce uncertainty, lower costs, and optimize production. Petroleum engineers, production engineers, geophysicists and geologists now use computerized models of the earth's crust to plan exploration and manage production of petroleum reservoirs.
Reservoir conditions are not static and gas and liquids in a reservoir can move rapidly. Engineers use numerical simulations in modeling the movements of gas and liquid in hydrocarbon reservoirs. These simulations aid engineers in gaining a better understanding of a reservoir's structure and flow properties, as well as aiding in developing an optimal production strategy which maximizes recovery and achieves the best economic result. Modern reservoir modeling has rapidly grown in complexity and can involve the integration of numerous software packages and computational engines. However, a modeling task can generally be broken down into five major steps: geological modeling, gridding and property generation, initialization, history matching and prediction.
Geological modeling concerns establishing reservoir geometry, reservoir boundaries and faults, as well as establishing basic rock properties such as porosity and permeability distributions. In order to provide an accurate representation of a reservoir, data can be incorporated from a very large number of sources including physical surveys, seismographic surveys, and data from wells. Engineers often have very detailed data about the characteristics and properties of specific and discrete locations within the reservoir. During initialization, this data may be used to determine the total amount of fluids as well as fluid compositions and phase saturations in discrete cells. History matching may be performed after running a simulation by comparing numerical results (e.g., for well production rate, gas oil ratio, water-cut, etc.) against actual field data and adjustments may be made to the reservoir model to reduce discrepancies. Once a model is refined, it may be used to predict future well production rates and a well management strategy based on the simulated reservoir flow conditions may be designed that aims to optimize overall recovery or economic results.
Gridding is typically performed to discretize a reservoir into a finite number of cells in order to simulate the behavior of the reservoir. It is often the case that grid orientation, grid size and grid geometry has a large impact on the accuracy of a simulation. In general, to resolve pressure, saturation, or permeability variation across the reservoir, smaller grid cells need to be used where such variation is large. For example, cell sizes in the surroundings of a well may need to be small because of the usually large pressure drawdown at those locations.
Besides grid sizes, grid orientation and grid geometry may also have a strong impact on simulation results. To reduce grid orientation effect, for example, grids may be aligned along major reservoir geological features, such as facies or faults, as well as streamlines for fluid flows. With commonly used Cartesian grid using local refinements, it may be possible to generate grids with cells of varying sizes, but it may be difficult to align cells along geological features or flow streamlines. For isotropic permeability distributions and features that are horizontal or vertical, this grid orientation problem with Cartesian grid may be partially solved with the advent of a K-orthogonal grids where gridnodes can be distributed more or less freely in space to conform to reservoir geometry. For anisotropic and highly heterogeneous permeability fields, it is almost impossible to create K-orthogonal grids. Therefore, conventional two point flux approximation (TPFA), which computes the flux of a face between two adjacent gridblocks by piecewise linearly interpolating the grid cell pressures, may no longer be accurate or valid and multi-point flux approximation (MPFA) should be used. MPFA involves a larger stencil to calculate the flux and is recognized by researchers as a way to improve numerical accuracy.
With MPFA, formulas for calculating fluxes are first derived for single phase flows, and then generalized to multiphase phase flows by accounting for saturation effects on relative permeability. The goal in the derivation is to solve for a pressure gradient field in the reservoir which yields fluxes that are consistent across interfaces between grid cells. To compute the pressure gradient in 2D, interaction regions are constructed around each vertex of the grid by joining the center of each grid cell sharing the vertex with centers of the cell edges that meet at the vertex. In 3D, interaction regions are created by constructing surfaces that are bounded by lines connecting block centers and face centers and lines connecting face centers and edge centers of cells around the vertex.
Computing a pressure gradient using MPFA may be described with reference to FIG. 1, which shows an exemplary set of adjacent grid cells 1101-1104 sharing a common vertex O. An interaction region 120 is constructed by joining cell centers C1-C4 with centers M1-M4 of cell edges that meet at the vertex O. To facilitate discussion, each portion of an interaction region that is contained in one of the cells is referred to as a subvolume. For example, in FIG. 1, quadrilaterals C2M1OM2, C3M2OM3, C4M3OM4 and C1M4OM1 are the four subvolumes of the interaction region 120.
With conventional MPFA, the pressure gradient is assumed constant in each subvolume and is obtained by using first order Taylor expansion or linear interpolation based on the pressure value at the cell center (Ci, i=1 to 4) and those at the edge centers (Mi, i=1 to 4) of the same subvolume. Using the pressure gradient and Darcy's law,F=nTλK∇p  (1)fluxes may then be written in terms of the pci and pek, where F is flux in feet per day, n is the normal vector, λ is the fluid mobility in 1/cp, K is a permeability tensor in millidarcys, ∇ is the gradient operator in 1/ft and p is pressure in pounds per square inch. Though not written explicitly, a unit conversion factor of 0.00633 may be used on the right hand side to convert md-psi/ft-cp to ft/day. By balancing fluxes across each half-edge inside the interaction region (M1M2, M2M3, M3M4, and M4M1), a formula for pressure at any point pek may be derived as linear functions of pci. Upon substitution, these equations may yield expressions for fluxes as linear functions of all pressure values at grid cell centers. The exact forms of these linear relationships may depend on the selected cell geometry and permeability tensors on grid cells surrounding the subvolume. In reservoir simulation, pressures at block centers may then be solved to determine the flows and state of the system.
In general, with MPFA the increase in computational cost when compared to TPFA is relatively small but results are shown to be more accurate, as expected. However, it has been observed that MPFA could result in severe unphysical oscillations in numerical solutions for single phase pressure equation−∇·K∇p=g  (2)if the permeability field is strongly anisotropic, where g is a source/sink term and a constant mobility of 1 is assumed to simplify the notations. These numerical oscillations may be linked to the lack of monotonicity in the solution matrix for the discretization problem and efforts have been made to improve MPFA so that matrix monotonicity can be restored. In any case, the possible occurrence of these oscillations may pose a severe limitation to the use of MPFA in many applications.
Accordingly, what is needed is an improved method for approximating flux gradients.